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INTilODUCTION 

Methods  of  solving  a  differential  equation  by   digital  conputer  are 
neither  ne\-r   nor  few.   The  question  of  the  neans  of  obtaining  the  solution 
icith  a  i.ia:'CLiaun  accuracy,  a  rrLniiiiui^  nuraber  of  computer  instructions  and 
storage,  and  a  iriiniauii  amount  of  calculation  tiiae  consistent  vxith  stability 
of  the  answer  has  been   the  source  of  e:>rtensive  investigations.   Standard 
numerical  analysis  texts  give  raany  procedures  but  these  are  rjst  the  subject 
of  this  report. 

This  report  reviews  pre^'/ious  works  on  trapezoidal  convolution  and 
Romberg  integration  for  digital  solution  of  differential  equations. 
Numerical  examples  of  solving  first  order  differential  equations  vjith 
above  methods  are  presented.  Efforts  have  been  made  to  iinprove  the  accuracy 
of  the  solution.  Simpson's  rule  and  higher  order  Newton-Cotes  quadratures 
require  ultimately  recursive  calculations.  Only  instantly  recursive  calcu- 
lations are  amenable  to  Z-transform  notation. 

It  is  necessary  to  note  that  trapesoidal  convolution  cannot  be  carried 
out  readily  on  a  digital  computer  because  of  the  large  number  of  multipliers, 
adders,  and  storage  capacity  required  if  the  operation  is  carried  out  in  a 
parallel  scheme  or  the  increasingly  long  computation  times  vxhen  done  serially'-, 
rnis  difficulty  is  by-passed  by   using  recurrence  relations. 


TRAPEZOIDAL  COWOLUTION 

Trapezoidal  convolution  can  be  used  to  solve  time -invariant,  tine- 
variant,  continuous  nonlinear,  discontinuous  nonlinear  equations  i-rith 
arbitary  initial  conditions.  There  are  many  ways  of  utilizing  trapezoidal 
convolution  to  solve  a  given  differential  equation.   Each  of  these  is  called 
a  program.  Three  classes  of  programs  are  discernible:   the  multiple-integrator 
substitution  program;  the  many-single-integrator  prograra;  and  the  single- 
integrator  program,  Halijak  (^)  has  generated  a  multiple-integrator 
substitution  program  and  a  single-integrator  prograra,  the  recursive  program. 

Approximate  Convolution 

The  solution  of  a  linear  differential  equation  t'Tith  constant  coefficients 
has  the  form  of  a  convolution.  The  digital  computer  must  approximate  convo- 
lution on  a  denumerable  set  of  equally  spaced  points.  The  Riemann  sum 
appro^djaation  is  the  one  commonly  used  on  digital  differential  analyzers 
and  finite  differences.  This  first  member  of  the  Newton-Cotes  quadrature 
family  is  by-passed  for  the  trapezoidal  quadrature. 

The  trapezoidal  approximation  proceeds  as  follows: 
rt  n-1   /"(i^+l  )T 

J    f(r)  g  (t-r)  dr=  X  /  ^^"^  s  (nT-r)  d^-        (i) 

0  k=o  At 

^1     Y.    i"(^<T)   g  (nT-kT)  +  f(kT+T)  g  (nT-kT-T)1=  J 
2     k=0  '^  ^ 

The  Z-transf orm  of  the  sequence  i  J    I  is 

X   J  2-  =  T(Zf)   (Zg)  -  0.5T|  f(0)Zg  +  g(0)Zf  =  Z(fg)  (2) 


n=0 


The  transforn  for  the  Simpson's  rule  approxLmation  need  not  be  derived 
because  it  and  higher  order  quadratures  do  not  yield  an  algebra  under  the 
operations  of  addition  and  discrete  convolution.   The  abbreviation  of 
trapezoidal  convolution  to  T-convolution  vri.ll  be  enployed. 

The  progran  for  solving  a  given  n-^h   order  differential  equation  is 
generated  by  viev^ing  it  as  the  integral  equation. 

y(t)  -r/^  g(t-r)  yWdT=  / ^  il=2ll.  :.:(r)dr        (3) 

^0  <)   (n-1)! 

where  g(t-T)  is  a  poljmonaal  of  degree  (n-1  )  and  xCt"/  is  the  forcing  function. 
It  is  evident  that  the  primary  task  is  to  find  the  approximate  Z-transforms 
of  l/s^  and  Y/s''^  in  order  to  solve  this  problem  on  a  digital  computer.  This 
can  be  done  appropriately  employing  equation  (2)  and  the  additional  starting 
fact  that  Z(l/s)  =  l/(l-z). 

Recursive  Prograai 

The  recursive  prograiTi  uses  a  sequence  of  ascending  order  differential 
equations  derived  from  the  given  differential  equation  to  solve  same. 
Consider  the  differential  equation 

(D^  H-aD^  -r  bB  +  c)  y(t)  =  ::(t),   D  =  d/dt, 

7(0)  =  y(0)  =  f(0)  =  0,   0  ^  t  ^  00  (^) 

First,  set  up  the  differential  equation 

(i>i-a)  y^(t)  =  1,   y^(0)  =  0  (5) 

The  Laplace  transform  of  this  D.E,  yields  after  division  by  s 

(l.|)y,  =1^ 


Talcing  the  Z- transform  of  both  sides  and  employing  the  approximate  Z-transform 
equation  (2),  of  the  convolution  ^  y,  yields 

S.I 

■  2y,    ^ 2Tz ^  V   f-J \  (6) 


= 2Tz ^  2:   j^-J \ 

[l-z]  [(2+aT)   -   (2-aT)zJ  ^s(s+a)) 


Tii3  above  appro:>CLriiation  is  now  used  in  the  next  derivate  y   (t),   vrhich  satisfies 

(D^-faT+b)  y^(t)  =1,       y   (O)  =  y  (O)  =  0  (?) 

2  2  2 

Dividing  now  by  s(s+a)  and  solving  yields 
Zy  i  Td+z)  2Tz 


2(1-z)  (2-i-aT)  -  (^-2bT2)2  +  (2-aT)z2  ^3^ 

=  zf 1 

(s(s^  +  as  +  b). 

The  third  derivative  yields  the  complete  solution  of 
ilP   -f  aD^  +  bD  +  c)  y(t)  =  x(t) 

y(o)  =y(o)  =y(o)  =  0  (9) 

Dividing  its  Laplace  transform  by  s(s  t-  as  +  b)  and  using  the  previous 
approximation  yields 


Zy  i t3  zd-i-z)  rzx-o.5x(o)3 


(10) 


(2+aT)  -  (6+aT-2bT^-cT^)z  +  (6-aT-2bT^+cT^)z2  -  {2-3.1) z^ 
Initial  conditions  other  than  zero  can  be  introduced  at  the  last  step. 
Tne  approximating  recurrence  relation  is  prepared  from  the  transformation: 
^  '  Zy  ~>  y         z^'Zy  ->  y^__^^  (11) 


mtegratoi-  Substitution  Program 


The  integrator  substitution  program  casts  the  Laplace  transfora  of  the 
given  differential  equation  into  inverse  powers  of  s  by  a  suitable  division 
and  substitutes  for  then  definite  functions  of  z.  First  of  all,  the  approxinate 
sampled  values  of  t^/n! ,  wliich  corresponds  to  the  transform  l/s   ,  are  four-d 

for  a  large  positive  integer  n.  The  derivation  of  the  approydinate  sampled 

r  / 
values  depends  on  the  fact  that  the  initial  value  of  t  "fnl   is  zero  except 

for  n=0.  For  non-urdty  positive  integers  n,  one  calculates  that 


The  above  is  now  used  as  a  descending  recurrence  relation 


s  s--i    I  s^-'v 


T(1+z) 
2(1-z)J 


'^ib^) 


1 \  \^h+zx 


s"  -/  UCl-z) 


(12) 


(13) 


until  it  stops  at 


^2(l-z)J   '  '^'■~^  ,  2,  3,  .  .  .  (1^) 


Calculating  Z(l/s~)  fror.i  equation  (2)  and  replacing  n  by  (n-l),  yield 


2/U  =  I2 .  r  T(l-t-z) 


n-2 


'(s^)   (l-z)2  I   2(i-z)J    .  ^--2.  3.  ^>,  .  .  .  (15) 

Of  course,  the  case  of  the  missing  integer  has  the  knovrn  value 

■  \l]  "  TTI  (16) 

which  cannot  be  continued  from  the  previous  formula.  One  should  note  tl-.at 
for  n=1  ,  2,  3»  "the  appro:djnate  Z-transfom  of  l/s'^  coincides  '-Jith  the  e:-:act 
Z-transfom.  The  integrator  substitution  program  is  at  hand  when  one  finds 
the  sampled  values  of  (l/s^)  f  (s).  '.vlien  n=l  ,  direct  calculation  i-rith  (2)  yields 


„  n  -\  ■  Td+z) ,-    Tf(o) 

VJhen  11=2,    3»    .    .    .»   direct  calculation  XiJith  equations   (I5)  and  (2)  yields 


lE(l-z), 


Differentiator  Substitution  Prograra 

Tiie   differentiator  substitution  program  does  not  require  an  intermediate 
Laplace  transfonaation,  but  does  require  one  more  initial  condition  which  is 
readily  obtained  from  the  given  initial  conditions  and  differential  equation. 

The  derivation  of  the  now  program  depends  on  equatiors  (I6)  and  (1?)  and 
the  Laplace  transforms  for  D  f  (t).  One  must  deal  trith  tx:o  categories  of 
derivatives,  the  first  and  the  others. 

The  approximate  Z- transform  of  the  first  derivative  is  now  obtained, 
5irst,  take  the  Laplace  transform  to  obtain 


Df(t)  =   sfCs)  -  f(0)  =  g(s)  (19) 

Divide  by  s  to  obtain 

li  =  f-lf(0)  (20) 

Sampling,  T-convolution  and  identification  of  g(0)  =  f(0)  yields 

T(±t4  2"  .  MQO  ^  2f  -  ^^  (21  ) 

zCUz)   ^^       2(1-z)   ^^   1-z  '^'^'  '' 

Solving  for  Zg  yields  the  desired 

^V^.xl  z.  2(l-z)  ^   2f(0)  -  T  f(0)  /ppN 

Z.Lx(o)i  -  ^rr-^   Zx  -   Vr(i,-2)  ^22; 

The  ne^rc.  category  of  derivatives  contains  those  of  order  two  and  liigher. 
Starting  ;^th 


^^^  =  S-^f(s)   -  Yl   ^"^'^  ^^"'^(0)  =  E(s).   i^^  2  (23) 


and  dividing  by  s^  yields 

Tald.n2  the  appi-o:>ci2i2.te  Z-transform  jd.elds 


Tz  (y  [zi  -  0.5  s(o)]  =  -^  -  z:  ^  (^ttt)'  ^^^^^ 


(25) 


Keeping  in  riind  the  transition  from  l/s  to  l/s     and  noting  that  g(0)  = 
f^'^'ho)  -lelds 

-  1  n-k-1      fv)  (n) 

-   1:  r^C'-^n  i"^    no)  +  0.5  f       (0).      n^  2  (26) 

::=1    L?C1-i-z)J 

The  presence  of  the  teras  f        (O)  in  equation  (22)   and  f       (0)   in  equation 

(26)  evidence  the  need  for  one  nore  initial  condition — a  very  slight  penalty 

for  elird-nating  the  inten.iediato  Laplace  transfon-iation. 

Tliis  differentiator  substitution  prograir.  is  equivalent  to  the  integrator 

substitution  and  recursive  prograia. 

IM??.0^/3D  T?.A.PSZOIDAL  II-ITSGPJITIOM  F0R:-IUIA 
A:^  its  SFF2CT  01!  ZIAI>3Z01DXL  COin^OLUTIO:i 

The  inprovenient  eiTiploys  an  alternative  fom  of  the  Euler-Xaclaurin 
fonr.ula.     ilecall  that  the  Taylor  series  representation  of  a  function  vjith 
reiiainder  is 

S(x)  =  g(0)  +  xg^^\o)  -i-  4  o^~^C0)  -i-  4  2^''*(0)  -^ /""  ^^^=^  g^''^(u)du 

(27) 
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If  the  given  function  is  integrated,  then  lengthy  computations  involving 
repeated  use  of  this  Taylor  series  yield 

1  r^  g(x)dx  =  ^^^  -  K  !0 til ^^^(1,)  ^28) 

h^Q  2     IE       2 


R(h)  =  )Lr    (A-  2>.^  H-A?')  g^^^  Uh)dX 


(29) 


The  remainder  Il(h)  is  given  by 

Since  the  polynomial  in  the  integrand  is  always  positive  on  the  interval  of 
integration,  the  first  mean  value  theorem  for  estimating  integrals  can  be 
applied  and  yields 

R(h):^  a—  g^^^  (9h)     0<  S<  1  (30) 

120 

The  essential  point  of  this  forraiila  is  that  the  added  terms  involve  only  the 
even  derivatives  of  the  function  at  the  ^rA  points  of  the  interval. 

Improved  trapezoidal  integration  utilizes  the  first  two  terms  of  the 
approximate  integral.  A  lengthy  application  of  this  approximation  to  the 
convolution  of  functions  f (t)  and  g(t)  yields 

Z(fg)  =  T  [(Zf)(zi)  -  0.5(fQ  zi  -  g^  Zf)j 
[(zf)(zi-)  -  o.5(fQ  zi+  g^  Zf)] 
^[(z^)(zl)  -o.5(fQ  zi  +  g^  Zf)) 
-  %  [(zf)(Zg)  -  o.5(f^  zg  H-  s^  zf)]  •  (31 ) 

If  this  formula  were  used  to  solve  a  differential  equation,  the  solution  of 

a  recurrence-differential  equation  vrould  be  required.  This  is  a  more  difficiilt 

problem,  than  the  original.  The  impasse  can  be  broken  by  using  the  trapezoidal 


12 


V 


dei'ivativG  substitution  prograai.  Hovrever,  difficulties  can  be  anticipated 
because  the  second  derivative  introduces  points  outside  the  integration 
interval. 

The  improved  trapezoidal  int32ration  in  conjunction  with  the  trapezoidal 
derivatives  leads  to  the  general  bilinear  forn 


Z(fi)  =  c^(Zf)(Zi)  -:-  ^Zf  -r  ^Z5  +  ^ 

Computation  yields 

°    ''  z(l-^z)2  J 


(32) 


(33) 


13  =  Tg.  p  -  26z  ~  122^  -  ^kJ  +  z'y  ^^^n  -  2z  +  3z^ 
'  "-  2^1-2   (i-i-22)  J  I  12   (i+z)2     - 

"  ^  2i^2   (l-z2)  -•  L      12   (1+z)2   J 

^  -  T      r  -1    +  7z  -^  2~  4-  z^l   ^  ^         T^  f  5  -  2z  H-  z^- 1  .     « 
12    I  z(l4-z)2  J   ^0"0   -2^        ^,^,)2      J  ^0^0 


(3^) 
35) 


-  22  "^  2" 
(1-^)2 


'P-^      r         ,7  "1     *     *  171J     •• 

t-  Lr^wJ  Vo^fe  Vo 


J  "0^=0 


(36) 


Tito   stumbling  blocks  are  apparent  in  these  formulas:      the  presence  of  the 

2 

factors  l/z  and  l/(l-j-z)  'trhich  represent  a  predictor  and  unbounded  function 

oscillating  about  zero  respectively. 

Recall  that  solution  of  a  differential  equation  reo^uires  knowledge  of  the 

sampled  values  of  l/s^  and  (l/s^)i.   The  first  step  is  to  find  Z(l/s'^).  It 

n+l        n+2  •>^— T 

is  possible  to  find  Z(l/s"  ),  Z(l/s*  ),  and  Z(l/s**'-^)  which  implies  that 

descending  recurrence  relations  can  be  found  which  decrement  the  e:cponent  by 


10 


one,  t'.7o,  or  three  integers.  The  long  series  of  calciilations  are  not 
e:-±ibited.  Hovrever,  the  recurrence  relation  which  decreraents  by  two 
integers  is  the  only  one  which  yields  the  expected  e:cact  Z(l/s  )  and 
Z(l/s  ).  This  determines  which  systeai  should  be  used  in  the  next  step. 

Calculations  of  ZCg/s*  )  produce  or.l.7  tvro   formulas  free  of  l/(l-i-z) 
factors.   These  are 


7  M  -\  ^  T  f1+sl  r  -1  +  I i!.z  -  2^  1  -   T(1  -  122  -  z^ 
^  is  ^J"  2  ll-zJl     iTi     J  2S  -   242  (1-z) 


'0 


4- 


T~(l-i-z)  „     ^'P 


2i^(1-z)  ''O   48  ^0  (37) 


7  f  1   "^  •  T^  f1  +  1Qz  4-.. 2^1  -  _  T^  1  +  5z  „  4-  T^   .    /^^N 


The  1/(1+2)  factor  or  its  powers  introduce  spurious  approiximations  even  for 

small  sain.pling  interval  sizes. 

_  J. 

Testing  the  formula  for  Z(g/s)  in  the  differential  equation  of  e 

cuic-ily  reveals  that  this  foriuula  is  spurious,  though  for  a  different  reason. 
The  predictor  l/z  is  the  source  of  difficulty. 

The  formula  for  Z(g/s  )  is  satisfactory  and  is  similar  to  the  Boxer- 
Thaler.  Z-form  operator  for  second  order  integration.  A  test  computation  in 
the  differential  equation  for  sin  a)  t  reveals  that  the  resulting  errors  are 
Sjialler  than  those  obtained  using  the  unimproved  integrator  substitution 
program.  On  the  other  hand,  some  test  computations  on  a  second  order 
differential  eq.uation  vrith  non  2ero  coefficient  for  the  first  derivative  term 
has  revealed  that  T-convolution  yields  better  approxLination  than  those  of 
Z-form.  It  is  therefore  conjectured  that  equation  (38)  is  best  enployed  in 
differential  eo^uations  involving  only  even  order  derivatives. 
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Use  of  an  ii'iprovsd  trapezoidal  quadrature  introduces  a  maze  of 
possibilities  vxliich  \Than  traversed  yields  a  lone  double  integration  formula  1 

R0:-S2RC-  i:JTEGIL'i.TIOrT 

Romberg  Algoritl^jn 

Consider  the  IJowton-Cotes  formulas  of  integration, 
(1  )     Zie  trapezoidal  approiomaticn  of  the  second  order 

/  f(x)dx  =  1  (f +f )  -  h  f^^hv*      B  =  1  (39) 

(2)     Simpson's  rule 

r1  H  .     B, 


'0 


f  (x)d.x  =  j  U^f.j^^  )-^T^  f^-"(5).         -'^  =  5J       (W) 


(3)     Gauss  formula 
•1 


J     fC:c)dx  =  ^  (7f^  ^  32f    .   -12f    ,  +32f       +7f  )   -  3    77  f^°'(S) 
-t  90         0  l/i^         1/2         3/^       1  2-^     6) 


3     =  ^  .      (^1 ) 


6       i^2 


vrhere  f  _  =  f  (x). 


One  can  observe  some  disadvantages: 

1  )  The  formation  of  tlie  coefficients  is  complicated,  and  it  is  not  e3.sj 
to  obtain  a  great  nuir.ber  of  these  coefficients  nor  to  calculate  them  in  a 
recurrent  marjrier; 

2)  Ey  increasing  the  number  01   points  for  a  fixed  interval,  one  needs  to 
calculate  again  other  numerical  values; 

3)  Some  of  the  coefficients  become  negative  for  n=3.  For  n=9  they  are 
all  positive,  but  for  n  :^10  there  are  negative  coefficients.  This  tends  to 
produce  ooor  roundoff  nro-Derties. 
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As  a  result,   the  higher  order  -Ie"'.-rbon-Cotes  foriaulas  are  seldon  used. 
The  order  of  the  error  ter^is   jurcps  by  2  in  going  frora  an  odd  nuraber  to  the 
ne:rt  even  nuraber,   which  tends  to  favor  the  even-order  fomulas. 

V/e  desire  an  algorithm  wliich  '-.'e  could  apply  easily  to  the  digital 
coviputcr;    thJ.s  algorithrr.  i;iust  be  such  that  the  coefficients  do  not  appear 
e:vplicitly;    in  other  words  >   it  vrill  not  ccnpel  us  to  calculate  the  numerical 
value  of  the  integrand.      Roriberg's  method   (10)  has  this  advantage  and  besides 
the  error  ter.^s  are  net  worse  than  those  of  the  ilewton-Cotes  formulas.     The 
coefficients  are  always  positive. 

To  obtain  his  algorithir-,   Roi'iborg  observed  that  Simpson's  rule  is  a 
linear  cor.ioination  of  trapezoiaaj.  formula:  tj 

1  (f     -i-  i^f    ,     ^-r    )   =±\l(b.^r  tlb]   -1-1    (Zlll^tl) 
6   ^^0  1/2       'V       3  U^2  2  ^       2    ^      2  2 /J 

"312         2  J 

i-Ioreover,  the  Gauss  formula  is  a  linear  coiiibination  of  the  Simpson's  rule: 

50^?^0^^^^lA-'^2^l/2-^32f3/^^7f^) 

15  12  I  5     6     6    )    -  2[     6     '       6      6/, 
.  t  f-0..^^-l/2  ;  -11 

"5  IT  "^^  Tj 


•7ne  combination  is  such  that  the  error  terms  disappear  and  the  order  of  error 
is  two  degrees  higher  than  that  of  the  Hewton-Cotes  o^uadrature.   Consequently, 

starting  iTLth.   trapezoidal  sums 

(2)   -       -^ 


n 


1  -^ 
?  "0 


.  1  -O 
-/n  2  lj 

n  =  1,  2,  ^,  8.  .  .  .  (i^3) 
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Ronber^  proposes  the  follovdns  algorithn: 


^n      =  2n     -^  Jl^ J} 

2n 
^   -  1 

v:here  ra  =  1  ,  2,  3i  -^-i  .  <>  . 

n  =  1,  2,A^  ,  .  .  .  i^) 

(2)    (2)         f2) 
lie   can  thsrofore  express  T   in  function  of  T    ,  T     ...  T^^  '   using 

2:n  1      2  '^*' 

2^  +  1  points,  ?or  ra  =  0,  1  ,  2  this  algorithn  produces  fornulas  which  are 
identical  vrith  trapezoidal  for::iula,  Sinpson's  r^ole,  and  Gauss  formula 
respectively,  but  for  n  5j  3  the  resulting  f orr^ulas  are  different  from  the 
Newton-Cotes  formulas. 
If  we  define 

(2)  =  1  ^  .p^^^    n  =  1.  2,  ^,  3,  .  ,  ,        (ij.5) 
the  evaluation  of  tho  trapezoidal  sums  becomes 


,  '2)  ^i  ^  (2)  _,.^   (2))     n  =  l.  2,4,  ,  .  .  •        (i^) 
2n     2   ^^      n 


and  the  higher  order  terms  vrill  be 

(2m-i-2)  _  ^^     (2m)       t     (2:.i)       -  (2m) 

^     2n  "   "r- 


2^^  2^ 


E^^''  -  1 


n^  m  (^7) 


1i^ 


The  whole  Romberg  algorithm  can  be  shox-m  bj  the  folloi-ring  diagran: 

(2)  ^(2)  .^(2) 

1  2  .  k 

^   (2)^     .   ^2)__^^   (2)__   \   (2) 

12  ^'  ^ 

\  (6;         \  (6) 


Diagonal  elements  of  high  order  are  not  necessarily  more  accurate  than  those 
of  lower  order;  this  has  been  first  shown  by  Bauer  (1 )  by  means  of  test 
calculations. 

In  the  same  manner  Romberg's  second  algorithra  defines  the  recurrence 
formula. 


,   (2m-f2)   .  (2m)   U  ^^^^^  -  U  ^^^"^ 


m  =  1 ,  2,  3,  .  .  . 

n  =  1,  2,1;.  ,  .  .  .  (48) 

If  ire  elaborate  the  diagram  of  T  and  U,  we  are  able  to  use  the 
following  relations 

2  n    =  -^ 


^  (2m+2)   ^2^m-1_^-^  ^  (2m)  ^  2^-ai-\    y  (2m) 

2        _  n n 

2^"  .  1  «50) 


m—  1,2,3»»»»       n  —  1,2,^-,.». 
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(2m+2)  (2)  (2) 

.'e  can  also  eroress  T  as  a  linear  coraoination  of  ?  ,    U 


^n  1      1 


(2)         ^2) 
U     .  .  .  u  ''   iTith  positive  coefficients. 
2         2a        ^ 


Nuraerical  Example 


To  illustrate  the  application  of  Romberg's  integration  process,  consider 


a  definite  integral 


I  =  /  f  cos  |:  ::  dx  =/  f  (:.:)dx 


^0         "       "0 

-•rhich  has  e^i^ct  solution  of  1.  ./o  first  divide  the  whole  interval  0,  1   in 

eight  ea,ual  parts,  then  using  equation  ('v5)  to  calculate  U  ,  U  ,  U, 

1         2        4' 

(2) 
U/       =  f(l/2) 


U,^^^  =\  [f(l/^)  +f(3A)] 

uj^^  =  I  [i(l/8)  -f  f(3/3)  H-  f(5/3)  -r  f(7/3:] 


Frora.  equation  (43) 

t/^^  =|[f(0)   -Hfd)] 

Then  using  equation  (46)  yields 

,  (2)  .  i  (,  (2)  ,  ,  (2)) 

2  2  ^^1       1    ^ 

3  2  ^4      4   ^ 

(4)        (4)        (4)        C6)        (6) 
Finall^r  enoloying  equation  (4?)   to  generate  T  ,   T  ,   T  ,   T.         ,   T  , 

-  2  4  d  4  8 

and  T     successively. 
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Kimerical  results  are  as  xollo/7s: 

T     ^     =  0.73539715 
1 


(2) 
U  ^  ^  =  1. 11 0721 3 

(2) 
T  =  0.9^305920 

2 

(2) 

u  ^  ^  =  1.0261722        u  '""'  =  0,99793920 


T   ^    '^  =  I.OC 


TT     (^) 


i22798 


T  ^"^  =  0.93711570 

u  ^^^  =  1.006^5^5 


(2) 


0.99673510 


(8) 


T   ^    ^  =  1.00013^5 

4" 

U  ^^^  =  0.99988200 

Hr 

T  =  1.0000082 

8 


Ilote  that  T  is  nost  accurate. 

8 

(2) 
Error  ?  =  -0. 00321^^9 

8 

Error  T  ^   ''  =   0.0000082 
8 

Error  T     =  -0.0000002 

8 

(8) 
Error  T     =  -0.00000007 


(6)  _ 


(6) 


^ 


=  0.99999150 


1.0000081 


(6) 


(8) 


T    =  0.99999980   T^'  '  =  0.99999993 
8  o 


ie-7T0N-C0TES  QUADRATURES  FOR  DIGITAL  SOLUTION 
0?  DI?FSREITTL\L  EQUATIOK^S 


The  use  of  integral  equations  to  establish  e:d.stence  theorems  is  a 
starjaard  device  in  the  theory  of  differential  equations,  both  ordinary  and 
partial.  It  oires  its  efficiency  to  the  smootliing  properties  of  integration, 
as  contrasted  TJith  coarsening  properties  of  differentiation.  If  two  functions 
are  close  their  integrals  raust  be  close,  whereas  their  derivatives  may  be  far 


1? 


apart  and  may  not  even  exist.  It  follows  that  a  key  to  tlie  solution  of  a 
differential  equation  lies  in  the  conversion  to  the  proper  integral  equation. 

Trapezoidal  Rule 

Dahlquist  (3)  has  shovm  that  tho  trapezoidal  rule  has  the  sjr^llost 
truncation  error  ainong  all  linear  i.iultistep  methods  "i-rith  a  certain  stability 
property.   Tlius,  trapezoidal  quadrature's  modest  accuracy  perforiiiance  is 
compensated  for  by  better  stability. 

Consider  the  vector  first-order  differential  equation 

^  +  A;7  =  V  (51  ) 

ax, 

Integration  yields 

W(t)  +  A  /       w(r)  d-r  =  V-     -:-  /      v(7-)  dr  (52) 

liiployins  trapezoidal  rule  yields 


:^  =  w 

0  0 


(53) 


W,     +   hi    (-,i       +  17    )    =   IT       +    i    (v       +   V    )  (54-) 


2        0  1  0       2       0  1 


K  rr\ 


w^  -^  ~  (v7     +  2ir     +  w  )  =  \:     +  -  (v     +  2v     -r  v  )  {^S) 

2        2        0  1  2  0       2       0  1  2 

Talcing  differences  of  adjacent  pairs  of  equations  yields  the  desired 

recurrence  relation 


^j     -  ^.r  "   +  *^  i:^     +  -^-^  .  )  =  -  (v  4-  V  J 

n    n-1  ^   -n         n-1    2   ^^    '^'^"' 

or 

1  -  AT  -  .    . 

w  =  Z  ,,    ^£_>Vl7n-1^ 

^   W  AT  n-1      ^--^7 (35) 

2  2 
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w  =2i-:i        +  ^  n>1  (57) 

n     n-1 

I'lote  that  equation  (57)  is  identical  Tvlth  the  recurrence  relation  obtained  by 
applying  trapezoidal  convolution  to  equation  (51).  It  is  worthwhile  to  notice 
tliat  if  the  coefficients  ^  and  B  are  functions  of  sarapling  interval  T  only,  then 

2 


n-i 


^^^^''^^■^^^r  ^d^  (^a^-^ ,  .  .  -.  Pol 

which  is  2iuch  nore  convenient  to  use  than  equation  (5?)  because  we  can  obtain 

the  desired  value  directly  xrithout  going  through  step-by-step  calculation. 

vrnen  n  approaches  infinity,  vt     converges  to 

n 

liv.i  w  =  __L„  (59) 

provided  that  such  lir^iit  exists. 

The  quadrature  iiriproverient  process  due  to  equation  (28)  improves  the 

second-order  integration  forrri-ala  onl^'.     If  an  additional  term  such  as 
(.'.  (-:•)  +  „  (^))/2 

wore  added  to  the  quadrature,   then  thJLs  weald  improve  the  fourth-order 
integration  formula,   and  so  on.      The  contradictory  demand  that  these  improved 
formulas  presents — that  not  only  the  past  but  the  future  values  of  the  unknown 
solution  m.ust  be  used — effectively  blocks  iciprovement  of  first-order  integration 
forrrdlas.      It  is  natural  to  turn  to  the  higher  order  NexTton-Cotes  q.uadi'atures 
for  imtjrovement  of  thjis  formula. 
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.liio'ii^-.cr  possibility  is  to  ir.iprovo  Sinpson's  rule  07  adding  higher  order 
terir.s.     ?.easo>iins  similar  to  that  employed  to  derive  equation  (28)  Halijak(6^ 
shows  that  a  :.iora  general  forr.  oi"  Sinpson's  rule  is 

i_[    cU)dx  =  ::o — -1    -2    ..v  r^o :js :_i2_., p,,) 


J 


.^(h)  ^  -■>—  s^^^  ( e  h)  c  -  0  ^  2  (60) 

1512 

3nployins  only  the  first  tenr.  ^ivos  Sinpson's  quadrature  rule.  This  result 
i-rill  not  be  used  for  accuracy  inprovenent  of  solutions. 

Appl;;,Hr.2  Sirapson's  rule  to  the  vector  first-order  differential  equation 
(51  )  L.d-elds 

w  =  w 

0  0 

vr_  +  ~  {\7    +  \r  ;  =  w^  +  -  (v^  -^  V  ; 

1  2   0    1     0   2   0    1 

2  ^  '  0     1    2     0   3   0     1    2 

3  I2  '  0       r      3  '  1    2   3  J 

=  W   +^  (v   -i-  V  )  -r  -  (v   +  ^V   +  V  J 

0  i-   0   1    ;-   1    2   3  '^ 
I'ote  that  tT-ro  estinates  of  the  sar,ie  integral  are  being  generated.  Tne  first 
cstiriato  occ-ors  vrith  the  odd  nunbered  equations,  and  this  cstinate  r.ust  be 
started  vjith  a  trapezoidal  integration.  The  second  estiriate  occurs  vrith  the 
even  numbered  equations. 

Talcing  differences  of  adjacent  pairs  of  even  numbered  equations  and  odd 
nuribered  equations  yields  the  syste:.i  of  equations 
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\-T    =    W 
0  0 

A  fyi  rji 

12^010201 

!■:     -  T-r         +  ~  (ij-     +  ^w     ,    ^-  tt       ) 
n         n-2        3        n  n-i  n-z 

=  -  (v     +  iiy         +  V       )  for  n  ^  2  (61  ) 

3       n  n-1  n-2  • 

Note  that  a  drastic  change  of  s^'-stematics  has  occurred.  The  calculations  are 

rec-orsive  orJ.y  for  n  ^  2.    ' 

Defirition:  A  system  which  is  recursive  after  a  finite  nuiaber  of  values  1^^ill 

be  called  an  ultimately  recursive  calculating  systera. 

Definition:  A  systen  which  requires  only  one  value  to  become  recursive  x^ll 

be  called  an  instantly  rec-orsive  calculation  system.. 

It  is  these  that  are  amenable  to  Z-transform  notation.  The  most  conspicuous 

example  is  trapezoidal  convolution. 

One  slight  change  in  the  system  of  equations  can  improve  accuracy.   If  . 

w  can  be  calculated  by  an  independent  method  to  a  greater  degree  of  accuracy 

than  by  trapezoidal  quadrature,  then  over  all  accuracy  i-7ill  be  improved.  A 

readily  available  method  is  the  Taylor  series  expansion  about  the  time  origin. 

The  resulting  ultimately  recxirsive  system  is 

w     =  w 
0  0 

2  ■   3 

•       1  0  0       2i      0       3j      0 

AT 

IT     -  vr     ^  +  -^-Cw     T  -^w         -f-  \i     ^ ) 
n         n-2    •     -5      n  n-1  n-^ 

=  :7(v     -f  4v        +  V       )     for  n  ^  2  (62) 

3     n  n-1  n-2 
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Follo'.vrins  the  sane  procedure  we  can   derive  ultimately  recursive  systems  for 
higher  order  Newton-Cotes  quadratures. 

The  corresponding  integrator  substitution  progran  does  not  exist  if  the 
conputer  menoiy  is  liirited  and  coraputor  tirae  r.ust  be  short. 

Fxmctional  Iteration 

The  investigation  starts  x-Jith  eo^uation  (3).  Simpson's  rule  is  then 
applied  and  the  resulting  equations  become  an  instantly  recursive  system. 
VJe  shall  call  this  metliod  "the  functional  iteration  method"  which  xrill  be 
denoted  as  ST.  From  equation  (3) 


y(t)  -:- /  g(t-r)  y  ('>)  d«r  =  /  g(t-r)  xCt)  dn- 


u 


y(t)  =  -  /   g(t-r)  y»dT-:-  Z^"  s(t-7-)  x:t)  d? 
^0  -<) 

we  generate  trapezoidal  solution  of  above  equation  first.  Then  the  desired 

functional  iteration  solution  t-rill  be  obtained  Ij/  applying  Simpson's  rule  or 

higher  order  Ne*.-rton-Cotes  quadratures  to  the  right-hand  side  of  the  folloiri.ng 

e  aquation: 

'^(t)  =  -J     g(t-r}  y:^.dr  -7  "  g(t->:  x('^;d7-      (63) 

where  y\t)   is,  as  mentioned  above,  the  trapezoidal  solution. 

Si-IPIRICAL  2RR0R  ANALYSIS 

A  comparison  of  the  computational  errors  generated  hy  various  sampling 
sizes  vrhen  using  T-convolution,  functional  iteration,  and  hi[;her-order 
Hewton-Cotes  quadratures  is  presented  for  the  follovring  first-order  ordir-ary 
differen':.ial  equations, 
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4r +  y  =  1  yCo)  =  o 


||-fy  =  0  y(0)   =   1 

^  +  ty  =  t  y(0)  =  2 

#  +  y^  =  1  y(o)  =  0 

In  order  to  carry  out  error  analysis  the  folloi-ring  noraenclatures  are 

employed: 

JE  ^{   ~   the  largest  absolute  value  of  error  in  the  fixed  period  of 

computation  (0.3  seconds); 

T  =  sasipling  interval  size;      . 

T  =  trar>ezoidal  convolution; 
c 

ST  =  trapezoidal  convolution  followed  by  Simpson's  rule; 
S  =  Sirapson's  l/3  2*ule; 
3/8  =  Simpson's  3/8  rule; 
G  =  G-auss'  quadrature. 

Derivation  of  .Recursive  Formulas 

The  process  of  derivation  of  the  rec^orsive  formula  for  the  differential 
equation  ^^-^7  =  0,  y(0)  =  0,  -.-rill  be  shovrn  in  detail.  For  the  other  three 
cases,  orHy   recursive  formulas  "i-jill  be  presented. 

1.  dv 

dt 

a)  T 


§?  H-  y  =  1       y(0)  =0  -(6^1^) 

do 


c 

1 
^7  "^  7  =  s 

-        1  -   1 
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Z'r  +  I 


(iFF)  =  My 


-y  ^  -    -~  -y  =  "'^ 


9 


2       \-Z  (l-2)~ 


2  1  -a 


(1    H.  |)Zy  -   (1    -  |).Zy  =  ^ 
J1    -^|)   -y^_^    (1    -|)=J0.    T.   T.   T.    .    .    .   } 


y 


y      = y         +0=0 

0  14-1         -1 


1 

- 

^ 

fji 

■~ 

1 

+ 

T 

^n- 

-1 

4-    ,  , 

" 

1 

T 

n 

2 

2 

y. 

— 

2. 

0 

- 

■T1 

m 

y 

<* 

4.  , 

2 

n   ;>    1  (63) 

Ilote  that  if  -^^e  apply  the  trapezoid  rule  alone  to  equation  (6^),    the 
sane  recurrenoe  forriula  as   {i^S)  "i-n.!!  be  obtained. 


-T  "  ~ 


t  =  ^oV  '^^-y  "^^  ^^°^ 


Y    =  y 
-^0      "^0 

y     =y     +T-|(y     +y)  (6?) 

10  2       0  1 

y     =  y     +  2T  -  I  (y     +  2y     +  y   )  (63) 

2  0  2       0  12 

Taking  the  difference  of  the  adjacent  formulas  we  obtain  the  recurrence 
foritiiila 


2i^ 


■^r   =  y   ^  +  'T"  _  £  (y   -r  y    ) 


1  -^ 

Tr                     -r- 

T* 

2 

n-1        1 

T     — 

2 

7     = 


9  _  T      ,   ?■; 


11   2  +  T  "^n-l  ■  2  -^  T  (69) 

which,  is  identical  vrith  equation  {6S) »     It  tiill  be  emphasized  that  i-fne-n 
er.iplo^'i.n^  the  trapezoidal  rule  alone  it  is  much  easier  to  get  recurrence 
formula  than  using  T-convolution  method  as  far  as  the  first-order  ordinary 
differential  equation  is  concerned,  but  the  same  is  not  true  for  higher- 
order  ones.  In  this  particular  case,  the  coefficients  of  y's  are  the 
function  of  ?  only,  so  equations  (53)  and  (59)  ^.re  applicable. 

>  n,     _   ./  n 


-n=^%^Mrf^J  =  P(f-^j 


(70) 


v'  -  2  -  T     r  -      2T 


'.-here   C(  =  -^^-^-i     /^  = 

2-^T,    r   2-i-T 

2T 

lira  y  =  -i —  =  HtlL.  =  1  (71  ) 

n-^co  n   1  -od    .,  _  2  -  T 


2  +  T 


b)  ST 


^2=^0'"''"3  ^^O-''^^  ^^-^ 


1    2 

->7i,  =  y„  -  ^T  -  i  (y  +  i^y     -f-  2y  -:-  ^,7  +  7  ) 
1^0        3   0     1     2     3    ^ 

?alclng  the  difference  of  the  above  eo^uations,  we  get  the  recurrence  formula 

>7  =   y?     ^  -i-  2T  -  ^  (y  -:-  4y    +  y   )  (72) 

/  -^         /n-2       3  ^-^n    "^n-l   '  n-2 

n  =  2,  ^>,   6,  3,  .  .  . 
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T-iiere  y's  are  trapezoidal  solutions, 
c)     3 

n 

V       =   -/       -i-    9'?    ^    -■    (^^       ~   ij-y       +   y     ) 


Tr  =   y       +   :j,;;:    _    -q-    (y       +   li-Y       -r    2y       +   -^7       "l-  Y     ) 
^^         "^0  '^  0  '^l  '^Z  '^3  "^I^ 

y.  =  7      +  2T  -  ?  (y  -^  'V  ,  ^y  .) 

i^         n-1  ^       n  -.-1  n-2 


J' 


r      -      6?         ^   3    _    ?   ^  i;/^ 

^-       3  t  ?    '   3  +  T  ^"-2  "  3T~T  -^n-l  n  ^   2  (73) 


I'ote  "ihat  instead  of  using  trapezoid  foriuula,   Taylor  series  eZ'^pansion  is 

en"oloyed  to  estimate  the  first  iterative  -ooint,  y  ,    to  ir^orove  accuracy, 

1 

d)     3/8 

1  0  Q       2i      0       JT     0 

..     =  y     +  (2T)  V     -^  (2?)"  y     .     Lgi:L.y  -l-  .    .    . 
.2  0  "0         21  °  3'1 

7    =7    ,  -^  3T  -  3/3T  (y    -i-  3y    .   -  3y    ^  +  y    J 
n         n-3  -^  n-i  n-2         n-3 

-^       3  +  3T       3  +  3?     -^^-3       3  -r  3T    *^n-1    '  •^n-2^ 

n:^  3  (7^) 

In  this  case  it  is  oh\a.ous  that  ::e  need  two  points  to   start  tlio  recursive 
forir.-ala  (7^)»     '"^e  shall  use  the  Taylor  series  to  goncrate  necessary  initial 
joints . 


\ 
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ii         n-^-  ^5         n  n-1  n-2  n-3  n-'-J- 


■^r      = 


130^ 


bd^JJ£l.r 


ST 


(8y 


3y    .  +  Sy       ) 
n-3 


n  ^    ^ 


(?5) 


Si:  +  y  =  0         7(0)  =  1 

Co 


o    _    f? 


n       2 


—  y 


^    1 


(76) 


.)     3: 


^n 


7n-. 


i   (7     +  ify  -f  7         ) 

3     "^n         '^n-l        ^^11-2 


n  =  2,   ^J-,   o,   b,    . 


(77) 


c) 


=    2^-^J:.   V 


i^T 


7 


3  +  T     ^"2       3  -H  T  "  ^-'' 


(73) 


d)     3/3 


7 


=  3^:^ 


n 


8  +  3T     n-3       8  +  3T     '  n-1        '  n-2 


(7         -^  7       )  ^  ^    3 


(79) 


y.    =  thi-nJii: 


fn    "^ 


'■O 


+  lifT     n-^       ^^5  -f  1 


T (87 


■'•■T 


n-i 


■^37         +  S7       ) 
^n-2  n-3 


n  ^    if- 


(30) 


0.7 

dt 


t7  =  t  7(0)  =  2 


2^nT2       ^n-1        ^  ^__  ^^^^ 


n    ^  1 


(31) 
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b)     ST 


n  =  2,  ^^,   6,  8,    .    .    .  (32) 


c)  S 

3  +  n^-  3  +  nT~  3  -  nT^     '^   '> 

d)  3/S 


^         8  +  3nT'-  '^  -^  ^"■■ 


y. 


n-3 


e) 


^"''''  f"    r.--1  W  -r    (r-2)Y  1 

3  +  3ni  ^ 


n  ^  3 


■v-'X 


!••  2Z   -H   y-    =    1  y(0)    =    0 

Q  o 

a)      - 


y.. 


-1    +       1    -  2T   (f-  y*^        -  7         -  T)  n  ^    1 


(8^) 


(35) 


(36) 


o;     S' 


n  =  2,   s'-,   6,   3,    .    .    . 


(37) 


23 


c;  i 


n  —  ~ 


^     3    ^3  n-1     n-2     n-?     '' 


;ar2. 


2T 
3" 


n  ^  2 


(83) 


d)  3/3 


\  = 


-1   v  /1   -  ^ 


■^  (37^^  ^  -  3y^ 
2  I  8     n-1      n-2 


+  y 


lz 


)  -  7   -3T 

n-3   . 


21 


n  S:  3 


(39) 


e;  G 


-1 


^v,  = 


51  (32y 


I2y 


■n-2 


-^  32y^    -^7-/       J  -J 


.-'■'^] 


28T 
^5 


n  5=  4 


(90) 


Results 


The  error  neasures  for  various  appro>d.inations  are  presented  in  Tables  1  , 
2,  3»  ^  i^-d  these  sane  error  neasures  are  graphed  in  Figures  1  through  11. 

It  should  be  pointed  out  that  the  values  of  error  raeasxires  are  presented 
vjith  floating  point  systen.  For  instance, 
0.112  i  -05  =  0.112  X  10"^ 
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Table  1.   j  1^^_,.| 


:or  p  -i-  y  -   I  ,  7(0)  =  0, 


.-- ( 

■:-coir7 

ij,  _' 

s 

->!  ^ 

:- 

0.003 

.1122-05 

.9302-06 

.3302-06 

.6502-06 

.1802-06 

0.00^:- 

.5602-06 

.5102-06 

.3902-06 

.4702-06 

.2602-06 

0.005 

.3302-06 

.1302-06 

.U:-02-06 

.4402-06 

.2602-06 

0.006 

.5SO2-O6 

.1602-06 

.2702-06 

.1802-06 

. 1 702-06 

0.007 

.9302-06 

.3702-06 

.3102-06 

.3502-06 

.2302-06 

0.003 

.1-^22-05 

.6902-06 

.1502-06 

.2302-06 

.1302-06 

0.009 

.2162-05 

.IO^a-05 

.1402-06 

.1402-06 

.1602-06 

0.010 

.2522-05 

.1302-05 

.1202-06 

.1602-06 

.1602-06 

0.012 

.3392-05 

.19^:'2-05 

.  1^^-02-06 

.1902-06 

.1302-06 

0.01i> 

.5692-05 

.2902-05 

.3002-07 

.1502-06 

.1902-06 

0.016 

.7332-05 

.3952-05 

.1802-06 

.1902-06 

.1102-06 

0.018 

.9':'72-05 

.^^922-05 

. 1 502-06 

.1402-06 

.1402-06 

0.020 

.1132-0^ 

.6302-05 

.1502-06 

.1602-06 

.8102-06 

0.022 

.1-^32-Oi^ 

.7-^32-05 

. 1 502-06 

.1602-06 

.9002-07 

0.02^1- 

.1702-0^ 

.8922-05 

. 1 302-06 

.1502-06 

.1202-06 

0.026 

.1992-0^ 

.1022-0^ 

.1372-06 

.1102-06 

.1302-06 

0.023 

.2322-Oii- 

.1202-C--> 

.1502-06 

.1302-06 

.1002-06 

0.030 

.2672-0^^ 

.1372-C^ 

.1  if  02-06 

.1502-06 

. 1 302-06 

0.032 

.3052-0^ 

.l6l2-0-> 

.1102-06 

.1102-06 

.1102-06 

0.03^ 

.3^32-0-^ 

.1762-04 

.1102-06 

.1002-06 

.1402-06 

0.036 

.3362-0^' 

.2022-0^ 

.1102-06 

.1002-06 

.1202-06 

0.033 

.^^312-0^ 

.2272-0^' 

. 1 082-06 

.1002-06 

.1002-06 

O.Oi^O 

^772-Oii 

.2552-C^ 

.1002-06 

.1002-06 

.9002-07 

0,0^2 

.5272-Oi!- 

.2772-C-. 

.1302-06 

.1002-06 

.1102-06 

O.OVi- 

.5772-0^ 

.3032-0^ 

.9002-07 

.1002-06 

.9002-07 

0.0^6 

.6302-0^ 

.3222-0^ 

.1302-06 

.1002-06 

.1502-06 

O.O-^I'S 

,63^^-2-0^:- 

.3^:42-0^ 

.1602-06 

.9002-07 

.1102-06 

0.050 

7^72-Oif 

.3932-0^ 

.6002-07 

.7002-07 

.1002-06 

0.052 

.8052-0^ 

.^llE-0^ 

. 1 202-06 

.9002-07 

.1002-06 

0.05^ 

.8622-0-^ 

.i^262-0^ 

.6002-07 

.1102-06 

.1002-06 

0.056 

93-:'2-0A'' 

.iu3ij,2-0^ 

.3002-07 

.6002-07 

.1302-06 

0.G53 

.9932-0^!- 

.ii'S52-04 

.3002-07 

.3002-07 

.3002-07 

Vj  •  W  s^  tj                            1 

.1072-03 

.5^62-0^ 

.lii.02-06 

.4002-07 

.1002-06 

0.062 

,1132-03 

.5^r32-0^ 

.1002-06 

.3002-07 

.1102-06 

0.06^ 

,1212-03 

.6132-oii' 

.3002-07 

. 1 1 02-06 

.1102-06 

0.066 

,1322-03 

.6332-0^ 

.3002-07 

.8002-07 

.9002-07 

0.063 

,1362-03 

.6552-0^^ 

.1002-06 

.5002-07 

.9002-07 

0.070 

.1^52-03 

.7272-0^^ 

.7002-07 

.9002-07 

.9002-07 

0.072 

,15^'^2-03 

.3022-0^ 

.9002-07 

.1002-06 

.1002-06 

0.07^ 

.1612-03 

.77^>2-0^- 

.1202-06 

.1702-06 

.3002-07 

0.076 

.1712-03 

.3502-0-> 

.1202-06 

.1302-06 

.9002-07 

0.073 

.1312-03 

.9322-0^ 

.1502-06 

.1402-06 

.6002-07 

0.030 

.1912-03 

.1012-03 

.6002-07 

.I5OE-O6 

.5002-07 

0.032 

.1972-03 

.9232-0^;- 

.1202-06 

.1402-06 

.9002-07 

0.03-.- 

.2052-03 

.1012-03 

.9002-07 

.1802-06 

.1102-06 

0.036 

.2202-03 

.1092-03 

.1002-06 

.2302-06 

.3002-07 

0.033 

.2312-03 

.1192-03 

.1502-06 

.2602-06 

.5002-07 
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Tabic  1  .      (cont.) 


l' 

T-conv 

oT 

i> 

3/8 

G 

0.090 

.236E-O3 

.1093-03 

.1303-06 

.2103-06 

.3003-07 

0.092 

.2^92-03 

.1133-03 

.1103-06 

.  2^1-03-06 

.1003-06 

0.09^' 

.2613-03 

.1233-03 

.1003-06 

.2703-06 

.1003-06 

0.096 

.2733-03 

.1373-03 

.  1 503-06 

.3203-06 

.9003-07 

0.093 

.2363-03 

.U^83-03 

.2203-06 

.3303-06 

,3003-07 

0.1  00 

.2993-03 

.1593-03 

. 1 603-06 

,^203-06 

.7003-07 

31 


blc  2. 


I  -..I 


7 


=  0,  7(0)  =1,    t  =  0.3 


_. . 

--conv 

ST 

--. 

_' 

:   I 

0.003 

.530"^-06 

.1013-05 

.6003-06 

.1003-06 

.1503-06 

. 1 503-06 

0.0C4 

.35o:i;-o6 

.9903-06 

.2603-06 

.1603-06 

.1^^03-06 

.2403-06 

0.005 

.135^-05 

.11 7-1-05 

.2203-06 

.1003-06 

.1503-06 

.2003-06 

C.006 

.1^3-05 

.1153-05 

.1^.^03-06 

.1303-06 

.1303-06 

.1303-06 

0.007 

.175--05 

.12^:'3-05 

.i6c:-j-06 

.1003-06 

.1002-06 

.3003-07 

0.003 

.2023-05 

.1^03-05 

. 1 ! 03-06 

.5003-07 

.1103-06 

.1703-06 

0.009 

c 2733-05 

.1733-05 

.  1  Ov'_j-0o 

.9003-07 

.7003-07 

.1203-06 

0.010 

.30i>3-05 

.1333-05 

.9003-07 

.9003-07 

.1003-06 

.1703-06 

0.012 

.;oo3-o5 

.2513-05 

.9003-07 

.5003-07 

.7003-07 

.1703-06 

.6133-C5 

.3^03-05 

.1503-06 

.6003-07 

.5003-07 

.2303-06 

0.016 

.7763-05 

.i^263-05 

.1103-06 

.3003-07 

.9003-07 

.3203-06 

0.01s 

.9793-05 

.5363-05 

. 1 1 03-06 

.1003-06 

.1003-06 

.4703-06 

0.020 

.t2i3-04 

.6593-05 

.1303-06 

.1003-06 

.6003-07 

.7003-06 

0.022 

.i^52-a'-> 

.7753-05 

.9003-07 

.8003-07 

.5003-07 

.3703-06 

0.02^' 

.1723-oi^ 

.9173-05 

.9003-07 

.6003-07 

■  .8003-07 

.1143-05 

0.026 

.2013-0^ 

.1053-Ci-;- 

.1103-06 

.7003-07 

.3003-07 

.1473-05 

0.023 

.23^-3-0^ 

.1232-Oii. 

.9003-07 

.;'-003-07 

.iK)03-07 

.1773-05 

0.030 

.2693-0^+ 

.1393-Oi^ 

.5003-07 

.7003-07 

.7003-07 

.2223-05 

0.032 

.3073-Oi^ 

.16^1-3-0.^ 

.1203-06 

.6003-07 

.5003-07 

.2653-05 

0.03^ 

.3^5:^-0^ 

.1733-0^ 

.3003-07 

.i:.003-07 

.i|O03-C7 

.3133-05 

0.036 

.3S73-0i> 

.20^3-0;'- 

. 1 003-06 

.6003-07 

.6003-07 

.3793-05 

0.033 

.^333-0i> 

.2293-0^ 

.3003-07 

.5003-07 

.5003-07 

.^4^>3-05 

OoOifO 

.4793-Oij- 

.2563-0^ 

.1203-06 

.^003-07 

.6003-07 

.5103-05 

0.0^2 

, 

. 1 003-06 

.5003-07 

.^003-07 

.5913-05 

0.0^.' 

.5792-0^ 

.30^3-0^!- 

.9003-07 

.^003-07 

.5003-07 

.6753-05 

O.Oi^ 

.6313-0-^ 

.32^3-0^ 

.1003-06 

.3003-07 

.5003-07 

.6393-05 

O.Oi^o 

.63^3-0i> 

'5''.<\->    oil 

.9003-07 

.6003-07 

.3003-07 

.3743-05 

0.050 

.7-:'9H;-0i^ 

.39o3-0v 

.5003-07 

.^■003-07 

.3003-07 

.4313-05 

0.052 

.8063-0^ 

.^>133-0i^ 

.1003-06 

.5003-07 

.5003-07 

.1113-04 

0.05^ 

. 

.1003-06 

.6003-07 

.6003-07 

.1243-04 

0.056 

.9363-0^ 

.i^'353-O;'' 

.7003-07 

.6003-07 

.^003-07 

0.053 

.9953-0^ 

.4363-0^.^ 

.3003-07 

.6003-07 

.5003-07 

.1533-04 

OoOcO 

.1073-03 

.5^73-0^ 

.3003-07 

.7003-07 

.6003-07 

.1703-04 

0.062 

.1133-03 

.5^93-Oii. 

.7003-07 

.9003-07 

.7003-07 

.1373-04 

O.Ooi^ 

.1213-03 

.61^3-0^ 

; 7003-07 

.6003-07 

.3003-07 

.2053-04 

0.066 

.1303-03 

.63i;-3-0i^ 

.9003-07 

.9003-07 

.2003-07 

.2243-04 

0.068 

.1363-03 

.6572-Oi^ 

.8003-07 

.1003-06 

.2003-07 

.2453-04 

0.070 

.1ii-53-03 

.7273-OiJ- 

.7003-07 

.1003-06 

.i+003-07 

.2663-04 

0.072 

.  155:^-03 

.8033-0^:- 

.1303-06 

.1203-06 

.3003-07 

.2903-04 

0.07'-^ 

.1613-03 

.7753-Oi^ 

.1303-06 

,1103-06 

.3003-07 

.3133-04 

0.076 

.1713-03 

0  "-■  •■•  ni. 

.6003-07 

.1503-06 

.i>003-07 

.3393-04 

0.073 

.1313-03 

.93i>3-0i^ 

.3003-07 

.1703-06 

.3003-07 

.3653-04 

0.030 

.1913-03 

. 1 01 3-03 

.6003-07 

.1303-06. 

.4-003-07 

.3933-04 

0.032 

.1973-03 

.9293-c^ 

.6003-07 

.2103-06 

.^^003-07 

.4233-04 

0.03i^ 

.2033-03 

.1013-03 

.S00£-07 

.2^03-06 

.4003-07 

.-•'43-04 

0.036 

.2203-03 

.1093-03 

.I303-O6 

.2303-06 

.5003-07 

.•  .73-04 

0.033 

.2313-03 

.1193-03 

.1^03-06 

.2803-06 

.4003-07 

.5203-04 

32 


■-:.-03 

o  o  r-  ""■'    ^  '^ 


/o. 


.6002-0? 
.cOOE^G? 


O533-0^ 

,^2M-Q^ 

.6312- 

-U'v 

.5?02. 

-0-> 

.7122. 

-C-^ 

7i^-"Vt;. 

...  1  z-^. 


>32CI.Qc 


_..^Cv  wO 


33 


jablG  3. 

W:l    '-  ^  -^- 

V  =  ^,  y(o)  =  2. 

0,    t  =2.0 

•■-1 

T-Gonv 

Cr.l 

. 

COOS 

.7IOE-O5 

. 1 932-04 

.2302-05 

.1302-05 

0.0C9 

. 5502-05 

.'  -^:-04 

.2002-05 

.1202-05 

0.010 

.^.!^-03-05 

J-Oi> 

.  ■;  ■^:-02-05 

.1102-05 

C.C12 

.3002-05 

:-05 

.1602-05 

.1002-05 

o.ou 

.6702-05 

.•.-..2-05 

.1302-05 

.1002-05 

O.Clo 

.1232-04 

.1-^:2-04 

.1402-05 

.6002-06 

O.OiS 

.1752-04 

]- 

.9002-06 

.8002-06 

0.020 

.2272-04 

•  -^  t' ,'-'  -'j^'^  ■* 

.9002-06 

<  -,  ^  -.      n^ 

0.022 

.2402-a4 

0    "n-l       r.  1, 

.1002-05 

»  _/  J  ^  _.  —  w    - 

C.02i^ 

.3522-04 

.4:':-2-0i> 

.6002-06 

.5002-06 

0.026 

.4192-04 

.5552-04 

.3002-06 

.5002-06 

0.023 

.4952-04 

.6612-04 

.7002-06 

.4002-06 

0.030 

.5732-04 

.7752-04 

.6002-06 

.4002-06 

0.032 

.6532-04 

.9002-04 

.9002-06 

.4002-06 

0.034 

.7502-04 

.1022-03 

.7002-06 

.4002-06 

0.036 

.3422-04 

.1152-03 

0.033 

.9402-0^!- 

.1302-03 

0.040 

.1042-03 

.1432-03 

0.042 

.1152-03 

.15:2-03 

0.04^' 

.1272-03 

.1762-03 

0.046 

.1392-03 

.1922-03 

0,043 

.1522-03 

.2092-03 

0.050 

.1642-03 

.2322-03 

0.052 

.I7S2-O3 

.2432-03 

0.054 

.1922-03 

.2712-03 

0.056 

.2072-03 

.2352-03 

0.053 

.2262-03 

.31 02-03 

0.060 

.2332-03 

.3322-03 

0.062 

.2552-03 

.3572-03 

0.064 

.2712-03 

.3792-03 

0.066 

.2892-03 

.4052-03 

0.063 

.3062-03 

.4262-03 

0.070 

.3252-03 

.4522-03 

0.072 

.3432-03 

.4702-03 

0.074 

.3632-03 

.5112-03 

0.076 

.3332-03 

.5332-03 

0.073 

.4042-03 

.55-2-03 

0.030 

.4252-03 

0.032 

.4462-03 

0.034 

.4692-03 

.6372-03 

0.036 

.4912-03 

.6342-03 

0.033 

.5142-03 

.7072-03 

0.090 

.5392-03 

.7532-03 

0.092 

.5622-03 

.7632-03 

0.094 

.31 52-03 

3^ 


Table  3-      Uont. ) 


0.096 
0.093 
0.100 

0.150 

0.200 


T-cor.v 


.6123-03 
.6393-03 
.oo5-3-03 

.1493-02 
.2663-02 


iT 


.0352-03 
.3903-03 
73-03 


r>,i  I  I 


,2023-02 
,3323-02 


3/8 
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Table  ^ 

ua.. . 

for  i^  +  r  =  n  7(0)  =  0,  t  =  0.3 

-1 

T-conv 

-i 

3/3 

- 

0.003 

.9212-03 

.2632-03 

0.00^ 

.1092-02 

.5352-03 

.3112-03 

./^732-03 

.^922-03 

0.005 

.^ii:-03 

.2192-03 

.iW72-C3 

.3^22-03 

.3162-03 

0.006 

.i!.50i^-03 

.1322-03 

.3672-03 

.2072-03 

.2652-03 

0.007 

.239:i:-03 

.1372-03 

.3532-03 

.1^52-03 

.1772-03 

0.003 

.3152-03 

.1212-03 

.1732-03 

.1232-03 

.1662-03 

0.009 

.13^3-03 

.8002-Oii. 

.lii-22-03 

.7332-0^ 

.1022-03 

0.010 

.5302-Oi^ 

.  1 032-03 

.1092-03 

.3iM2-0i^ 

0.012 

.1092-03 

.4732-0^^ 

.3072-0^1- 

.6272-0^ 

.5512-0^ 

0.01^ 

.6692-0^1- 

.2752-0^'4. 

.7732-0^ 

.3362-0^ 

.5l^S-0i{' 

0.016 

.5552-04. 

.2562-0^ 

.6122-0^ 

.2992-0^ 

..^762-oi;- 

0.013 

.6332-0^ 

.2372-04 

.3572-0^ 

.2472-0^ 

.3352-Oi^ 

0.020 

.5502-Oi^ 

.26l2-C-'4- 

.3922-0^!- 

.1672-0'^:- 

.2352-Oi^ 

0.022 

.■^5S2-0i> 

.2362-0-^ 

.3382-Oi; 

.2132-0^ 

.2152-0.4' 

0.02^ 

.6^52-0-'v 

.3122-0-!- 

.2132-0^ 

.1/^92-0^ 

.  1 5^2-0-^ 

0.026 

.'^532-0'^ 

.2252-0^ 

.  1 632-0-1^ 

.1002-0^ 

.2122-0^^ 

C.023 

.5952-0'^ 

.2332-0.^ 

.2l62-0i^ 

.1562-0^ 

.  1 1 22-0^1- 

0.030 

.6l32-0^> 

.2962-0-'> 

.1672-0^ 
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DISCUSSION 

Investigation  of  the  graphs  shows  that  there  e:d,sts  a  saiupling  interval 
sise  t-j-hich  iraniinizes  the  largest  absolute  value  of  error  for  the  fixed 
conputation  duration  in  T-convolution  (T  )  and  function  iteration  (ST). 
For  tho  case  of  ^  +  y   =  0,  y(0)  =  1,  the  optimura  sarapling  interval  size  are 
not  shoT-m  but  tlie  author  has  confirmed  this  by  other  coraputation  wliich  shows 
that  the  optim-un  interval  size  falls  between  0.002  and  0.003  seconds.  There 
is  sorae  he--u.-*istic  justification  for  the  existence  of  this  riiininuni  error  in 
the  T-convolution  case.  If  the  interval  size  is  large,  then  trapezoidal 
integration  generates  large  errors.  If  one  observes  the  function  of  z  obtained 

by  the  T-convolution  process,  then  the  liimit  as  T  approaches  zero  yield  a 

,    \n 
factor  C1-z;  in  the  denoriiinator,  where  n  is  the  order  of  the  differential 

equation  being  studied.  In  spite  of  the  exact  solution  being  bounded, 

unbo-oiided  approxLnating  functions  then  occur  when  r.  ^  2  and  the  errors  are 

correspondingly  large. 

The  mininuDi  error  point  is  a  function  of  coraputation  time  t.  Referring 
to  Figure  8,  one  i-rill  observe  that  for  T-convolution  \-ihen   computation  time, 
t,  increases  from  1  second  to  2  seconds,  the  minimum  error  point  moves  from 
T  =   0.01  to  T  =  0.12.  Few  examples  do  not  justify  any  definite  conclusion, 
but  from  the  above  example  the  rd-nimuai  error  point  may  be  the  function  of  the 
solution  time.  If  one  increases  the  computation  time,  t,  further,  one  can  find 
the  optimum  sampling  range  for  ininimum  error  point  for  the  particular  differential 
equation  under  consideration. 

Functional  iteration  technique  is  primarily  designed  for  improvement  of 
accuracy  in  computation.  Test  computations  show  that  it  surely  lowers  the 
maicLmuirL  error,  |E,„,.],  for  some  interval  of  computation  tirae  t,  but  as  t 
increases  further,  the  rate  of  increase  in  l^^^^l  exceeds  that  of  T-convolution 


^ 


and  it  "Till  becorie  less  accurato  than  that  of  T-convolution.  This  accuracy 
regression  property  of  functional  iteration  technio^ue  is  vividly  displayed  in 
Figures  7  sjnd   8.   In  this  particular  exanple,  when  T  lies  betiieen  0,012  and 
0.015,  the  accuracy  of  3T  is  greater  than  that  of  T  ,  but  for  _  >  0.016  or 
■_  "^  0.012  the  accuracy  is  poorer  than  T  .  Again,  there  is  a  problem  of 
estirtiating  optimum  sampling  range.  As  in  the  case  of  the  opti^ron  sampling 
ran^e  for  nininur.i  error,  one  can  only  assume  that  one  of  the  deciding  factors 
Tn.ll  be  computation  tiir.e,  t.  Since  the  accuracy  inprovenent  property  of  the 
functional  iteration  (ST)  is  not  hi^h  even  if  one  uses  the  sampling  size,  ", 
in  optir.ura  range,  unless  one  can  estimate  regression  point  beforehand,  one 
sho-uld  avoid  employing  this  method. 

Generally  speaking,  application  of  Simpson's  rule  and  other  higher-order 
"e'.rton-Cotcs  integration  formula  in  the  solution  of  differential  equation 
gives  better  results  than  ?-convolution  method.  Especially  \:hen   the  ir-itial 
values  are  estimated  as  accurately  as  possible,  the  solutions  arc  a::iazingly 
acc'^-ate.  The  effect  of  initial  values  to  the  accuracy  of  solution  is  clearly 
sho".."n  in  Figure  v.  '.7aen   trapezoidal  rule  is  used  to  approj±matc  y.  ,  the  error 
curve  of  Simpson's  rule  is  siirdlar  to  the  T-convolution,  thou:jh|E..,j,  |  of  ih.e 
for.::er  is  lovrer  than  that  of  the  latter.  On  the  other  hand,  tl.c  adoption  of 
Taylor  series  in  approxiiiiation  of  y  gives  a  much  better  error  characteristic. 
There  is  one  often  overlooked  trouble  that  occurs  when  Simpson's  formula  is  used 
to  start  the  chain  for  the  odd-numbered  sample  points.  I^ne   result  of  this 
jumping  t'.ro  steps  ahead  is  that  the  accumulated  errors  at  the  odd-  and  even- 
numbered  points  are  rather  independent  of  each  other,  especially  as  the  weights 
attached  to  the  computed  integrand  values  are  different.  This  tends  to  produce 
an  oscillation,  vrnile  this  oscillation  is  due  to  errors  committed  and  hence 
gives  some  measure  of  the  accuracy  of  the  results,  it  can  be  very  arjioying. 
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Oscillation  is  observed  in  everj   case  of  test  computations;  it  is  conjectiared 
that  for  special  case  of  y'  =  f(t),  oscillation  may  be  avoided  if  we  estimate 
y  hy  Tarflor   series ,  then  even-ntudbered  point  chain  only  is  computed  by 
Sir.ipson's  rule  and  the  half  Simpson's  formula  (Hamming,  p«  12?)  is  used  to 
produce  each  odd-numbered  point. 

Considering  that  Simpson's  second  rule  and  Gauss'  formula  have  to 
appro:>ajaate  more  starting  points,  that  coefficients  are  more  complex,  that 
acc-Liracy  improvement  effect  is  about  the  same  as  Simpson's  rule,  higher-order 
iTet'rton-Cotes  o^uadrature  formulas  have  doubtful  performance  in  digital  solution 
of  differential  equations. 
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Tiie  essential  content  of  tliis  report  is  laotaprograiXTiing  for  the  digital 
approiciiiiation  of  the  solution  of  first-order  difforontial  equations  using 
Newton-Cotes  quadratures. 

'.Jhen   trapezoidal  quadratures  are  employed,  there  are  two  possible 
programs:  the  first-order  vector  differential  eo^uation  prograra  and  the 
integrator  substitution  pro gran.  These  can  be  succinctly  described  in 
Z-transfonn  notation. 

Rom.berg  proposes  an  algoritlira  for  nucaerical  integration  which  generates 
high-order  integration  quadratures  by  repeated  use  of  linear  combinations  of 
trapezoidal  quadrature.   The  coefficients  of  Roraberg  o^uadratures  are  always 
positive  and  the  error  terns  are  not  worse  than  those  of  the  Ilewtoii-Cotes 
quadratures.  It  could  be  very  useful  in  digital  solution  of  differential 
equations. 

VJhen  Sinpson's  rule  or  higher-order  rIewton-Gotes  quadratures  are  employed, 
the  phenoraenon  of  ultimately  rec^orsiveness  appears  in  the  first-order  vector 
differential  equation  program,  but  there  is  no  integrator  substitution  program 
possible  for  high-speed  calculation.   Z-transform  notation  must  be  dropped  in 
ultimately  recursive  calculations. 

Empirical  error  analysis  shows  that  there  exists  an  optimuiu  sampling  . 

range  which  riiinirriizes  the  largest  absolute  value  of  error  for  the  fixed 

commutation  duration  in  T-convolution  (T  )  and  function  iteration  (ST);  that 

c 

as  a  result  of  accuracy  regression  property,  ST  has  only  limited  use  in  the 
solution  of  differential  equations;  that  Simpson's  rule  or  higher-order 
riev; ton-Cotes  quadratures  has  higher  accuracy  but  less  stable  than  T-convolution; 
that  the  estimation  of  starting  points  plays  a  vital  role  in  accuracy 
improvement  of  the  solution  of  differential  equations  by  higher-order  Nevrton- 
Cotes  cuadratures. 


In  conclusion,  if  one  employs  Nev7ton-CotGS  quadratures  to  solve  differential 
eo^uations,  one  can  expect  that  trapezoidal  rule  has  lo;:er  accuracy  but  higher 
stability  and  less  coiiputation  time  while  the  reverse  is  true  for  higher-order 
quadratures. 


